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Abstract 

Reconstruction of the CMB in the Galactic plane is extremely difficult due to the dominant foreground 
emissions such as Dust, Free- Free or Synchrotron. For cosmological studies, the standard approach 
consists in masking this area where the reconstruction is not good enough. This leads to difficulties 
for the statistical analysis of the CMB map, especially at very large scales (to study for e.g., the low 
quadrupole, ISW, axis of evil, etc). We investigate in this paper how well some inpainting techniques 
can recover the low-^ spherical harmonic coefficients. We introduce three new inpainting techniques 
based on three different kinds of priors: sparsity, energy and isotropy, and we compare them. We show 
that two of them, sparsity and energy priors, can lead to extremely high quality reconstruction, within 
1% of the cosmic variance for a mask with Fsky larger than 80%. 
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1. Introduction 

1.1. CMB and the masked data problem 

As component separation methods do not pro- 
vide good estimates of the Cosmic Microwave 
Background (CMB) in the Galactic plane and at 
locations of point sources, the standard approach 
for CMB map analysis is to consider that the data 
are not reliable in these areas, and to mask them. 
This leads to an incomplete coverage of the sky 
that has to be handled properly for further anal- 
ysis. This is especially true for analysis methods 
which operate in the spherical harmonic domain 
where localization is lost and full-sky coverage is 
assumed. For p ower spectrum est imation, methods 
like MASTER ( |Hivon et al.|2QQ2| ) solve a linear ill- 
posed inverse problem allowing to deconvolve the 
observed power spectrum of the masked map from 
the mask effect. 

For non-Gaussianity analysis, many approaches 
have been proposed to deal with this missing data 
problem. Methods to solve this problems are called 
gap filling or inpainting methods. 

The simplest approach is to set pixel values in 
the masked area to zero. This is sometimes claimed 
to be a good approach by assuming it does not add 
any information. However, this is not correct, as 
setting the masked area to zero actually adds (or 
removes) information, which are also not the true 



CMB values. A consequence of this gap filling tech- 
nique is the creation of a lot of artifactual coeffi- 
cients at the border of the mask in a wavelet anal- 
ysis or a leakage between different multipoles in a 
spherical harmonic analysis. This effect can be re- 
duced using an apodized mask. A slightly more so- 
phisticated method consists in replacing each miss- 
ing pixel by the average of its neighbors and iterat- 
ing until the gaps are filled. This technique is called 
diffuse inpainting a nd has been used in P lanck-LFI 
data pre-processing - - - 



Zacchei et al. 



( |2Q11^ . It IS ac- 
ceptable for treatment of point sources, but is far 
from being a reasonable solution for the Galactic 
plane inpainting in CMB non-Gaussianity analysis. 
In lAbrial et all (|2007|); lAbrial et all (|2008|), 



the problem was considered as an ill-posed mverse 
problem, z = Mx, where x is the unknown CMB, 
M is the masking operator, and z the masked 
CMB map. A sparsity prior in the spherical har- 
monic domain was used to regularize the problem. 
This sparsity-based inpainting approach has been 
successfully used for two different CMB studies, 
the C MB wea k lensing on Planck simulated data 
dPerotto et al 1|2Q1Q| [Plaszczynski et al.||2Q12[) , and 



the analysis of in tegrated Sachs- Wo lfe effect (ISW) 
on WMAP data (Dupe et al.J |2Qll[ ). In both cases, 
it was shown using Monte- Carlo simulations that 
the statistics derived from the inpainted maps can 
be trusted with high confidence level, and that 



sparsity-based inpainting can indeed provide an 
easy and effective solution to the large Galactic 
mask problem. 

It was also shown that sparse inpa inting is use- 
ful fo r weak lensing data (Pires et al.| 2QQ9), Fermi 
data dSchmitt et al. 2010), and asteroseismic data 
(Sato et al. 2010'). The sparse inpainting success 



has motivated the community to investigate more 
inpainting techniques, and other approaches have 
been recently proposed. 



Bucher fc Lomsl ( |2012[ ); |Kim et al.J ( |2012[ ) seek 
a solution which complies with the CMB proper- 
ties, i.e. to be a homogeneous Gaussian random 
field with a specific power spectrum. Nice results 
were derived, but the approach presents however 
the drawback that we need to assume a given cos- 
mology, which may not be wise for non-Gaussianity 
studies. For large scale CMB anomalies studies. 



Ben-David et al. (2012D, build 



de Oliveira-Costa Tegmark| 



1, 
Te 



Feeney et al.| ( |201lD ; 
ing upon the work o: 

(2006) proposed to use generalized least-squares 
(GLS), which coincide with the maximum likeli- 
hood estimator (MLE) under appropriate assump- 
tions such as Gaussianity. This method also re- 
quires to have the input power spectrum. In ad- 
dition, when the covariance in the GLS is singular, 
the authors propose to regularize it by adding a 
small perturbation term that must also to ensure 
posit ive-definiteness. Using this regularization, the 
estimator is no longer a GLS (nor a MLE). 

Which model for the CMB 7 

As the missing data problem is ill-posed, prior 
knowledge is needed to reduce the space of can- 
didate solutions. Methods attacking this problem 
in the literature assume priors, either explicitly or 
implicitly. If we put aside the zero-inpainting and 
the diffuse inpainting methods which are of little 
interest for the Galactic plane inpainting, the two 
other priors were: 

• Gaussianity: the CMB is assumed to be a ho- 
mogeneous and isotropic Gaussian random field, 
so it makes sense to use such a prior. In prac- 
tice, methods require the theoretical power spec- 
trum, and it has eith er to be estimated u sing a 
method like TOUSI ( |Paykari et al.||2012D , or a 
cosmology has to be assumed, which is even a 
stronger assumption than the Gaussianity prior. 
We should also keep in mind that the goal of 
non-Gaussianity studies is to check that the ob- 
served CMB does not deviate from Gaussianity. 
Therefore, we should be careful with this prior. 

• Sparsity: the sparsity prior on a signal consists 
in assuming that its coefficients in a given repre- 
sentation domain, when sorted in decreasing or- 
der of magnitude, exhibit a fast decay rate, typ- 
ically a polynomial decay with an exponent that 
depends on the regularity of the signal. For the 
CMB, its spherical harmonic coefficients show a 



decay in 0{£~'^) up to £ around 900 and then 
0{i~^) for larger multipoles. Thus, the sparsity 
prior is advocated, and this explains its success 
for CMB-related inverse problems such as in- 
paint ing or component separation ( ,Bobin et al. 



2012) 



1.2. Contributions 

In this paper, we revisit the Gaussianity and spar- 
sity priors, and introduce an additional one, namely 
the CMB isotropy, for the the recovery of the spher- 
ical harmonic coefficients at low £ (< 10) from 
masked data. We describe novel and fast algorithms 
to solve the optimization problems corresponding 
to each prior. These algorithms originate from the 
field of non-smooth optimization theory and they 
apply efficiently to large-scale data. We then show 
that some of these inpainting algorithms are very 
efficient to recover the spherical harmonic coeffi- 
cients for < 10 when using the sparsity or energy 
priors. We also study the quality the reconstruc- 
tion as a function of the sky coverage, and we show 
that a very good reconstruction quality, within 1% 
of the cosmic variance, can be reached for a mask 
with sky coverage larger than 80%. 



2. CMB Inpainting 

2.1. Problem statement 

Suppose that we observe z = Mx, where x is a 
real-valued centered and square-integrable field on 
the unit sphere S^, and M is a linear bounded 
masking operator. The goal is to recover x from 
z. 

The field x can be expanded as 



^=0 m= 



where . 



§2 



where the complex- valued functions Yirn are the co- 
called spherical harmonics, £ is the multipole mo- 
ment and m is the phase ranging from —£ to £. 
The a^^rn are the spherical harmonic coefficients of 
X. In the following we will denote S the spheri- 
cal harmonic transform operator and S* its adjoint 
(hence its inverse since spherical harmonics form 
an orthobasis of L2(S^)). 

If is a wide-sense stationary (i.e. homoge- 
neous) random field, the spherical harmonic coeffi- 
cients are uncorrelated. 

If moreover the process is isotropic, then 



2 



where is the angular power spectrum, which de- 
pends solely on i. 

In the rest of the paper, we consider a finite- 
dimensional setting, where the sphere is discretized. 
X can be rearranged in a column vector in R"^, and 
similarly for z G M^, with m < and a G C^. 
Therefore, M can be seen as a matrix taking val- 
ues in {0, 1}, i.e. M G ^mxn({0, 1}). The goal of 
inpainting is to recover x from z. 

2.2. A General Inpainting Framework 

The recovery of x from z when m < n amounts to 
solving an underdetermined system of linear equa- 
tions. Traditional mathematical reasoning, in fact 
the fundamental theorem of linear algebra, tells us 
not to attempt this: there are more unknowns than 
equations. However, if we have prior information 
about the underlying field, there are some rigor- 
ous resul ts showing that su ch inversion might be 
possible (IStarck e^ al.||2010). 

In general, we can write the inpainting problem 
as the following optimization program 

find X e Argmin R{x) s.t. z - Mx G C , (1) 

where is a proper lower-bounded penalty func- 
tion refiecting some prior on a?, and C is closed con- 
straint set expressing the fidelity term. Typically, 
in the noiseless case, C = {x : z = Mcc}. This is 
the situation we are going to focus on in the rest 
of the paper. As far as R is concerned, we will con- 
sider three types of priors, each corresponding to a 
specific choice of R. 

3. Sparsity Prior 

Sparsity-base d inpainting h a s bee n proposed for 



where | 



^llg = Z^ikiT^ where \z\ = 

for z e C. For q = 0, the Iq 
pseudo-norm counts the number of non-zero entries 
of its argument. The inpainted map is finally recon- 
structed diS X = S*a. 

Solving ([2| when q = is known to be NP-hard. 
This is further complicated by the presence of the 
non-smooth constraint term. Iterative Hard thresh- 
olding (IHT), described in Appendix A, attempts 
to solve this problem. It is also well-known that 
li norm is the tightest convex relaxation (in the £2 
ball) of the Iq penalty (Starck et al. 2010). This sug- 
gests solving (2| with g = 1. In this case, the prob- 
lem is well-posed: it has at least a minimizer and 
all minimizers are global. Furthermore, although 
it is completely non-smooth, it can be solved effi- 
ciently with a provably convergent algorithm be- 
longing to the family of proximal splitting scheme s 



longmg to tne lamny 01 proximal spnttmg scneme s 
( [Combettes fc Pesquet||2011| [Starck et aL]|2010| ). 
This can be done using the Douglas- Rachford (DR) 



algorithm described in Appendix B. 



Comparison between IHT and DR 

We have compared the IHT and DR sparsity-based 
inpainting algorithms on 100 Monte-Carlo simula- 
tions using a mask with sky coverage Fsky = 77%. 
In all our experiments, we have used 150 iterations 
for both iterative schemes, /3 = 1 and = 1 (Vn) 
in the DR scheme (see Appendix B). For each in- 
painted map i G {1, • ' ' 7 100} 7 we computed the 
relative mean squared-error (MSE) 



^True ^(0 



the C MB in Abrial et al. (2007); Abrial et al. and the its version in percent per i 



(2008) and for weak-lensmg ma ss m ap reconstruc- 



tion in 



Pires et al.| (2009 



( |2010| ); |Dupe et al.| poTT] ); [Rassat et al.| ( |2012p 71t 

was shown that sparsity-based inpainting does not 



2010). In Perotto et al. 



El 



100 X / e^'^ 



(%). 



sparsity-based mpamti: 
destroy CMB weak-lensing, ISW signals nor some 
large scale anomalies in the CMB, and is there- 
fore an elegant way to handle the masking problem. 
Note that the masking effect can be thought of as 
a loss of sparsity in the spherical harmonic domain 
because the information required to define the map 
has been spread across the spherical harmonic basis 
(leakage effect). 

Optimization problem 

If the spherical harmonic coefficients a of x (i.e. 
a = Sx) are assumed to be sparse, then, a 
well-known penalty to promote sparsity is the Iq 
(pseudo- or quasi-)norm, with q G [0, 1]. Therefore, 
becomes 

find a G Argmin ||a||^ s.t. z = MS*a , (2) 

aeCP 



Fig. [T] depicts the relative MSE in percent per 
i for the two sparsity-based inpainting algorithms 
IHT (/o) and DR (h). We see that at very low 
/i -sparsity inpainting as provided by the DR algo- 
rithm yields better results. We have performed the 
test with other masks, and arrived at similar con- 
clusions. As this paper focuses on low-^, only the 
li inpainting as solved by the DR algorithm will be 
considered in the sequel. 



4. Energy Prior 

Optimization problem 

If we know a priori the power-spectrum (C£,m)£,m 
(not the angular one, which implicitly assumes 
isotropy) of the Gaussian field x, then using a max- 
imum a posteriori (MAP) argument, the inpaint- 



3 



20 



Inpainting Error (Fsky = 77^ 



15 



^ 10 
o 



Sparsity (DR) 



■-■ IHT 




6 



10 



Multipole 



Figure 1. Relative MSE per I in percent for the two sparsity-based inpainting algorithms, IHT (dashed 
blue line) and DR (continuous blue line). 



ing problem amounts to minimizing a weighted ^2- 
norm subject to a linear constraint 



find X = argmin \\Sx\\ 



s.t. z Mx , (3) 



where for a complex- valued vector a, ||a||^_i = 
I 1^ 

m ^^"^ , i.e. a weighted £2 norm. By strong 
convexity, Problem (|3| is well-posed and has ex- 
actly one minimizer, therefore justifying equality 
instead of inclusion in the argument of the mini- 
mum in (|3|. 

Since the objective is different iable with a 
Lipschitz-continuous gradient, and the projector 
Proj|3,.^^]Vi3,| on the linear equality set is known 
in closed-form, one can use a projected gradient- 
scheme to solve (l3|. However, it turns that the 
estimate of the descent step-size can be rather 
crude, which may hinder the convergence speed. 
One can think of using an inexact line-search but 
this will make the algorithm unnecessarily compli- 
cated. This is why we propose in Appendix C a new 
algorithm, based on the Douglas-Rachford splitting 
method, which is easy to implement and efficient. 

In all our experiments in the rest of the paper, 
we have used 150 iterations for this algorithm with 



/3 = 50 and = 1 (Vn) (see Appendix C). This 
approach requires the power-spectrum as an 
input parameter. In practice, an estimate Q of the 
power spectrum from the data using MASTER was 
used, and we set C^^rn = The latter is reminis- 
cent of an isotropy assumption which is not neces- 
sarily true. 



5. Isotropy Prior 

5.1. Structural constraints on the power spectrum 

Strictly speaking, we would define the set of 
isotropic homogeneous random processes on the 
(discrete) sphere as the closed set 



where C^^rn = ^ 



\{Sx), 



and Ci is the angular 



power spectrum, which depends solely on £. Given a 
realization of a random field, we can then naively 
state the orthogonal projection of x onto Ciso by 
solving 



mm 



V - X 
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This formulation of the projection constraint set where, 
is not straightforward to deal with for at least the 
following reasons: (i) the isotropy constraint set in- 
volves the unknown true Q (through the expec- 
tation operator), which necessitates to resort to Proj<^ («£m) 
stochastic programming; (ii) the constraint set is ^ 
not convex. 



if \ai^rn\ > fJ^i^e 

if \ai^rn\ < IJ'i-S 

ai^rn otherwise. 



(5) 



5.2. Projection with a deterministic constraint 



An alternative to 
be to replace C^,^ 



the constraint set Ciso would 
n and by their empirical 



sample estimates, i.e. C^^rn by |a£,rn| and by 

= 2Wi^rn\^^.rri\^ • Howcvcr, this hard con- 
straint might be too strict in practice and we pro- 
pose to make it softer by taking into account vari- 
ability inherent to the sample estimates, as we ex- 
plain now. 

In a nutshell, the goal is to formulate a con- 
straint set, where for each the entries a^^rn of the 
spherical harmonic coefficient vector a deviate the 
least possible in magnitude (up to a certain accu- 
racy) from some pre-specified vector /l^; typically 



we take = Vw or its empirical estimate y C^. 
Put formally, this reads 

Ce = {x : a = Sic, ||a£,m| - ^ \/i,m} 

Cs is a compact set, although not convex. 



^9(a'€,m) 




3?(a^, 



Figure 2. The set 



We now turn to the projection on C^. We be- 
gin by noting that this set is separable, i.e. Cs = 

where is the band depicted Fig. [2] It turns 
out that for fixed /l^, is so-called prox-regular 
since its associated orthogonal projector is single- 
valued with a closed-form. Indeed, the projector 
onto is 



Projc^(a^) = S* (Proj^^ (a^,m))^ 



(4) 



Choice of the constraint radius: To devise a 
meaningful choice (from a statistical perspective) 
of the constraint radius e, we first need to derive 
the null distributior|^ of \ai^rn\ — M^7V(^, m), with 

the particular case when ii^ = \J~Ci. a^^ra being the 
spherical harmonic coefficients of a zero-mean sta- 
tionary and isotropic Gaussian process whose an- 
gular power spectrum is Q, it is easy to show that 



2 \ap, 



'C^X'(2) and 



2Lax,'^Qx'(2L), V^>2, 



where L = 2^ + 1. By the Fisher asymptotic for- 
mula, we then obtain 



3 1 
4' 4 



and 



1 - 



> 2 



where means convergence in distribution. Thus, 
ignoring the obvious dependency between |a^,rn| 
and /L^^, we get 



\Ci£,r, 



J_ L±l\ 
AL-> AL j 



Note that the distribution of this difference can be 
derived properly using the delta method in law, but 
computing the covariance matrix of the augmented 
vector (a£,m)^7 remains an issue. The quality of the 
above asymptotic approximation becomes better as 
I increases. 

The upper and lower critical thresholds at the 
double-sided significance level ^ a ^ 1 are 



e± = \/Ci 



where <l> is the standard normal cumulative distri- 
bution function. 

We depict in Fig |3] the upper and lower crit- 
ical thresholds normalized by y/C~i at the classi- 
cal significance level 0.05 as a function of I. One 
can observe that the thresholds are not symmet- 
ric. This entails that different values of e should 
be used in (^. Furthermore, as expected, the two 
thresholds decrease in magnitude as I increases. 



^ That is the distribution under the isotropy hypoth- 
esis. 



5 





Figure 3. The upper (left) and lower (right) normalized critical thresholds at the significance level 0.05 
as a function of I. 



They attain a plateau for I large enough (typi- 
cally > 100). It is easy to see that the two limit 
values are ^3/4 - 1 + 1.96/^8 = 0.559 and 
73/4 - 1 - 1.96/ V8 = -0.8269. 

The isotropy-constrained inpainting algorithm 
is given is Appendix D. 



6. Experiments 

6.1. ai^rn Reconstruction 

In this section, we compare the different inpaint- 
ing methods described previously: the DR-sparsity 
prior inpainting, the DR-energy prior inpainting 
and the isotropy prior inpainting. We also test 
the case of no-inpainting, applying just an Fsky 
correction to the a^^rn spherical harmonic coeffi- 
cients of the map (i.e. data a^^rn are corrected 
with a correction term equal to l/y^Fsky). We 
have run 100 CMB Monte-Carlo simulations with 
a resolution corresponding to nside equal to 32, 
which is precise enough for the large scales we 
are considering. We have considered six differ- 
ent masks, with sky coverage Fsky ranging in 
{0.98, 0.93, 0.87, 0.77, 0.67, 0.57}. 

Fig [4] displays these masks. Note that point 
source masks have not been considered here, since 
they should not affect a^^rn estimation at very low 
multipole. Each of these simulated CMB maps has 
been masked with the six masks. Fig. [5] top right 
shows one CMB realization masked with the 77% 
Fsky mask. Fig. [5] middle and bottom show respec- 
tively the inpainting of the Fig. [5] top right image, 
with the four methods (i.e. using sparsity, energy, 
isotropy priors, and Fsky correction). Fig. [6] depicts 
the results given by the ^l sparsity-based inpainting 
method for the six different masks. We can clearly 
see how the quality of the reconstruction degrades 
with decreasing sky coverage. 



The a^^rn coefficients of the six hundred maps 
(100 realizations masked with the six different 
masks) have then been estimated using the three 
inpainting methods and the Fsky correction meth- 
ods. The approaches were compared in terms of 
the relative MSE and relative MSE per see end 
of Section [3] for their definition. 

Fig. [7| shows relative MSE for ^ = 2 to 5 versus 
the sky coverage Fsky. The three horizontal lines 
correspond to 1%, 5% and 10% of the cosmic vari- 
ance (CV). Inpainting based on the sparsity and en- 
ergy priors give relatively close results, and anyway 
better than the one assuming the isotropy prior or 
the Fsky correction. It is interesting to notice that 
a very high quality reconstruction (with 1% of the 
CV) can be obtained with both the sparsity- and 
energy-based inpainting methods up to ^ = 4, for a 
mask with Fsky larger than 80%. If WMAP data do 
not allow us to make non-Gaussianity studies with 
such a small mask, Planck component separation 
will certainly be able to achieve the required qual- 
ity to make possible the use of such a small Galactic 
mask. With a 77% coverage Galactic mask, the er- 
ror increases by a factor 5 at ^ = 4 ! 

Fig. [8] shows the same errors, except that they 
are now plotted versus the multipoles for the six 
different masks. The ^1 sparsity-based inpainting 
method seems to be slightly better than the one 
based on the energy prior, especially when the sky 
coverage decreases. 

6.2. Anomalies in CMB maps 

One the main motivations of recovering all sky 
maps is to be able to study statistical properties 
such as Gaussianity and Statistical Isotropy on 
large scales. Statistical Isotropy is violated if there 
exists a preferred axis in the map. Mirror parity, i.e. 
parity with respect to reflections through a plane: 
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Mask (Fsky=98%) 



Mask (Fsky=93%) 





Figure 4. Six masks with respectively Fsky G {0.98,0.93,0.87,0.77,0.67,0.57}. The masked area is in 
blue. 



X = X — 2(x • n)n, where n is the normal vector to 
the plane, is an example o f a statistic where a pre- 



ferred axis can be sought ( [Land fc Magueijo||2005 
[Ben-David et al.||201 2). 



With all-sky data, one can estimate the iS-map 
for a given multipole in spherical harmonics, by 
considering ( Ben-David et al.]|2012 ): 



Seih) = ^ (-1 



(6) 



where a^^(n) corresponds to the value of the a^^n 
coefficients when the map is rotated to have n as 
the z-axis. Positive (negative) values of ^'^(n) cor- 
respond to even (odd) mirror parities in the n di- 



rection. The same statistic can also be considered 
summed over all low multipoles one wishes to con- 
sider (e.g. focussing only on low multipoles as in 
Ben-David et al.„2012J : 



S'tot(n) = ^ ^'^(n) 



(7) 



It is convenient to redefine the parity estimator 
as ^(n) = 5tot(n) - (Cax - 1), so that {S) = 0. 

The most even and odd mirror-parity directions 
for a given map can b e considered by estimating 
( |Ben-David et al.||2012[ ): 

max(5) - ij{S) 



S+ 



a{S) 



(8) 
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Simulated CMB (liiiax=10) 



Masked Simulated Data (Fsky=77%) 





DR Sparse Constraint Inpainting: Mask Fsky = 77% 



DR Energy Constraint Inpainting: Mask Fsky = 77% 






Isotropy Constraint Inpainting: Mask Fsky = 77% 



Fsky Correction: Mask Fsky = 77% 



§^ 




Figure 5. Top, input smoothed simulated CMB map (^max = 10) and the same map, but masked 
(Fsky = 77%) and not smoothed (i.e. input simulated data). Middle, inpainting of the top right image 
up to ^max = 10, using a sparsity prior (left) and an energy prior (right). Bottom, inpainting of the top 
right image up to ^max = 10, using an isotropy prior (left) and a simple Fsky correction (right). 



I min(6') — /i(5')| 



(9) 



where \i{S^ and cr(5') are the mean and standard 
deviation of the S map. The quantities 5+ and S- 
each correspond to different axes n+ and n_ re- 
spectively and are the quantities which we consider 
when discussing mirror parity in CMB maps. 

In order to use inpainted maps to study mir- 
ror parity, it is crucial that the inpainting method 
constitutes a bias-free reconstruction method, i.e. 
that it does not destroy existing mirror parities nor 
that it creates previously inexistent mirror parity 
anomalies. 



In order to test this, we consider three sets of 
1000 simulated CMB maps using a WMAP7 best 
fit cosmology with nside = 32, one set for each in- 
painting prior. For each simulation we calculate the 
mirror parity estimators S± . We isolate the anoma- 
lous maps, defined as those whose S± value is larger 
than twice the standard deviation given by the sim- 
ulations. This gives us 34 maps (3.4%) for the even 
mirror parity and 40 maps (4%) for the odd mirror 
parity in the full-sky simulations. 

Since recent work show s a potential odd- mirror 
anomaly in WMAP data ( |Ben-David et al.||2Q12| , 
we focus on testing for potential biases in oda- 
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DR Sparse Constraint Inpainting: Mask Fsky = 



DR Sparse Constraint Inpainting: Mask Fsky = 93% 





DR Sparse Constraint Inpainting: Mask Fsky = 87% 



DR Sparse Constraint Inpainting: Mask Fsky = 77% 






DR Sparse Constraint Inpainting: Mask Fsky = 67% 



DR Sparse Constraint Inpainting: Mask Fsky = 57% 




Figure 6. i\ sparsity-based inpainting of a simulated CMB map with different masks. Top, Fsky = 98% 
and 93%, middle 87% and 77%, and bottom 67% and 57%. 



mirror anomalies. Figj9] show the percentage dif- 
ference between the true S- value of the one esti- 
mated from the inpainted map. For this statistical 
test, the result is not dependent on the inpainting 
method as the three different priors give similar 
biases and error bars. For positive-mirror anoma- 
lies, we find similar results (not shown in Figure). 
Similarly to the previous experiment, we can see 
that a very good estimation of the parity statistic 
can be achieved for masks with Fsky larger 0.8. 

7. Software and Reproducible Research 

To support reproducible research, the developed 
IDL code will be released with the next version 



of ISAP (Interactive Sparse astronomical data 
Analysis Packages) via the web site: 

http : / / www . cosmostat . org 

All the experiments were performed using the 
default options: 

— Sparsity: 

Aim = cinb_lowl_alin_inpainting(CMBMap, 
Mask, /sparsity) 

- Energy: 

Pd = inrs_powspec( CMBMap, Mask) 
P = inrs_deconv_powspec( pd. Mask) 
Aim = cmb_lowl_alm_inpainting (CMBMap, 
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Figure 7. Relative MSE (in percent) per I versus sky coverage for different multipoles. The three hor- 
izontal lines correspond to 1%, 5% and 10% of the cosmic variance. Reconstruction within the 1% of 
the cosmic variance can be achieved up to ^ = 4 using the inpainting methods based on the sparsity or 
the energy priors for a mask with Fsky larger than 80%. With a 77% coverage Galactic mask, the error 
increases by a factor 5 at ^ = 4. 



/Energy, Prea=P) 

— Isotropy: 

Pd = inrs_powspec( CMBMap, Mask) 
P = inrs_deconv_powspec ( Pd, Mask) 
Aim = cinb_lowl_alin_ inpainting (CMBMap, 
/Isotropy, Prea=P) 

8. Conclusion 

We have investigated three priors to regularize 
the CMB inpainting problem: Gaussianity, spar- 
sity and isotropy. To solve the corresponding min- 
imization problems, we have also proposed fast 
novel algorithms, based for instance on proximal 
splitting methods for convex non-smooth optimiza- 
tion. We found that both Gaussianity and ^l spar- 
sity priors lead to very good results. The isotropy 
prior does not provide as good results. The spar- 
sity prior seems to lead to slightly better for > 2, 
and more robust when the sky coverage decreases. 
Furthermore, unlike the energy-prior based inpaint- 
ing, the sparsity-based one does not require a power 
spectrum as an input, which is a significant advan- 
tage. 

Then we evaluated the reconstruction quality 
as a function of the sky coverage, and we have seen 



that a high quality reconstruction, within 1% of 
the cosmic variance, can be reached for a mask 
with Fsky larger than 80%. We also studied mirror- 
parity anomalies, and found that mask with Fsky 
larger than 80% also permitted nearly bias-free re- 
constructions of the mirror parity scores. Such a 
mask seems unrealistic for WMAP data analysis, 
but it could and should be a target for Planck com- 
ponent separation methods. To know if this 80% 
mask is realistic for Planck, we will need more real- 
istic simulations including ISW effect, and residual 
foregrounds at amplitudes similar to those achieved 
by the actual best Planck component separation 
methods. Another interesting study will consist in 
comparing the sparse inpainting to other methods 
such as t he maximum likelihood or the Wiener fi l- 
tering ( |Feeney et al.||2011| |Ben-David et al.||2012 ). 
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Figure 8. Relative MSE in percent per i versus i for the six masks. 



Appendix A: Algorithm for the Iq problem Appendix B: Algorithm for the h problem 

Solving ([2| when g = is known to be NP- 
hard. This is further complicated by the presence 
of the non-smooth constraint term. Iterative Hard 
thresholding (IHT) attempts to solve this problem 
through the following scheme 



^n+l _ ^H^ (^n + s(z - MS*a^)) , (10) 

where the nonlinear operator is a term-by-term 
hard thresholding, i.e. A^(a) = (p^(ai))i, where 
p^{ai) = ai if |ai| > A and otherwise. The 
threshold parameter decreases with the itera- 
tion number and is supposed to play a role similar 
to the cooling parameter in simulated annealing, 
i.e. hopefully it allows to avoid stationary points 
and local minima. 



It is now well-known that li norm is the tightest 
convex relaxation ( in the £2 ball) of the Iq penalty 
dStarck et "aL]|2QlQ[ ). This suggests solving with 
g = 1. In this case, the problem is well-posed: it 
has at least a minimizer and all minimizers are 
global. Furthermore, although it is completely non- 
smooth, it can be solved efficiently with a provably 
convergent algorithm belonging to th e family of 
proximal splitting schem es (Combettes fc Pesquet 
"20111 [Starck et al.||2QlQD . 



In particular, we propose to use the Douglas- 
Rachford (DR) splitting scheme. Let /3 > 0, 
(<^n)nGN ^ sequcuce in ]0,2[ such that 

X]t^i^<^n(2 — an) = +00. The DR recursion, ap- 
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Figure 9. Percentage difference between anomalous (i.e. high) odd mirror parity scores {S-) before (i.e. 
fuh-sky) and after inpainting for three different inpainting priors. 



phed to ([2| with g = 1, reads 
a^+^ = Proj|^,,=MS*a} K), 

(11) 

prox^l|.|l^ is the proximity operator of the /i-norm 
which can be easily shown to be soft-thresholding 

prox^l|.|l^(a) = AS(a), (12) 

where A^(a) = (p^(ai))i and Px{ai) = 
sign(ai)max(0, la^l - P). Proj|^^^_MS*a} is 
the orthogonal projector on the corresponding 
set, and it consists in taking the inverse spherical 
harmonic transform of its argument, setting its 
pixel values to the observed ones at the cor- 
responding locations, and taking the forward 
spherical harmonic transform. It can be shown 
that the sequence {a^~^^)n G N converges to a 
global minimizer of ([2]) for g = 1. 



Appendix C: Algorithm for the I2 problem, 

This scheme is based again on Douglas-Rachford 
splitting applied to solve problem ([3|. Its steps are 

• =Proj|^^,^Ma.}(^"), 

-cc^) , (13) 

where P and an are defined as above, and the prox- 
imity operator of the squared weighted £2-norm is 

P^ox^l|S(.)f^_,(a^) =S* ((Sa^)® (^^)) ' 

(14) 

where stands for the entry-wise multiplica- 
tion between two vectors, and the division is 
also to be understood entry-wise. The projector 
Proj|3,.^^]y[3,| has a simple closed- form and consists 
in setting pixel values to the observed ones at the 
corresponding locations, and keeping the others in- 
tact. It can be shown that the sequence {x^~^^)n^fq 
converges to the unique global minimizer of ([3|. 
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Appendix D: Algorithm for Isotropy 
Inpainting 

The isotropy-constrained inpainting problem can 
be cast as the non- convex feasibihty problem 

find X eCsr\{x \ z = mx} . (15) 

The feasible set is nonempty since the constraint 
set is non-empty and closed; the CMB is in it under 



the isotropy hypothesis. To solve (15), we propose 
to use the von Neumann's methodof alternating 
projections onto the two constraint sets whose re- 
cursion can be written 

= Proj{,,,=M.} (Projc. (^("^)) • (16) 
Using closedness and prox-regularity of the con- 



straints, and by arguments from |Lewis fc Malick 



(^2008), we can conclude that our non-convex al- 
ternating projections algorithm for inpainting con- 
verges locally to a point of the intersection H : 
z = Mx} which is non-empty. 
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